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Abstract 

Suspensions of colloidal hard particles are a well established model system to study fundamental 
questions of statistical mechanics experimentally and theoretically. Their popularity as a model 
system is partly due to the assumption that solvent effects on phase transition kinetics can easily 
be factored out - at least if one considers processes that are very slow compared to diffusive time- 
scales. In particular, crystal nucleation rates for hard sphere suspensions are often compared in 
units of a "typical" diffusion time. 

Here, we present a study on the crystallization process in supersaturated suspensions of hard 
spheres, taking into account the solvent viscosity via hybrid dynamics simulations. We show that, 
although the induction times are long compared to the diffusion time, crystal nucleation rates still 
depend strongly on the solvent viscosity. We do not observe any differences in the sizes, shapes or 
structures of the nuclei that are formed. The effect is purely kinetic. As crystal nucleation kinetics 
are influenced non-trivially by hydrodynamic interactions, colloids might not be as ideal a model 
system to study phase transitions as has been assumed up to now. 
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Colloids are widely used as model systems to study fundamental questions of statistical 
mechanics. Over the past twenty-five years the phase behaviour, phase transition kinetics 
and glass transition of colloidal suspensions have been observed in numerous experiments 
and modelled by means of theory and simulations 1-7 . One reason the statistical mechanics 
community is interested in colloids, is that the interactions in colloidal suspensions can be 
tailored such that their equilibrium phase diagrams resemble those of atomic and molecular 
substances. Going beyond equilibrium physics, it is often suggested that also the phase 
transition behaviour on the colloidal scale should be similar to that on the atomic scale - at 
least in cases where excluded volume effects dominate the phase behaviour-"—. 

However, while atoms and colloids resemble each other in many aspects, there is a major 
difference in their dynamics: Colloidal particles are suspended in a solvent. They interact 
directly with each other (e. g. by excluded volume) as well as indirectly by means of mo- 
mentum transfer via the solvent. The latter phenomenon, called hydrodynamic interaction, 
is well known and thoroughly studied in the context of colloidal flow and sedimentation^. 
But when colloids are used as model systems for phase transition kinetics or for the glass 
transition, hydrodynamic interactions are usually neglected^^^. 

In this article we would like to address the effect of the solvent on crystal nucleation in 
suspensions of colloidal hard spheres. We chose this system, because there has been a long- 
standing debate on the large differences in crystal nucleation rates that have been observed 
experimentally as well as by means of computer simulation^-^^. 

Experimentally, hard sphere suspensions are synthesized in various ways. Common systems 
are polystyrene spheres suspended in water, and sterically stabilized polymethylmethacrylate 
(PMMA) spheres in an organic liquid such as decalin. When a suspension is prepared for a 
crystallization experiment the static properties of the solvent are carefully adjusted (i. e. the 
density and the refractive index of the solvent are matched to those the colloidal particles), 
and the concentration is determined as accurately as possible. Transport properties such as 
the solvent viscosity, however, are hardly ever mentioned in publications on crystal nucleation 
in colloidal suspensions. There seems to be an underlying assumption, that crystallization 
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kinetics of chemically different systems can be compared (irrespective of solvent viscosity), 
as long as one uses a reduced time-scale - reduced with respect to a "typical" diffusion time. 

There are several candidates for this time-scale. Some authors use the invers of D Q , the self- 
diffusion constant for an infinitely diluted system^^^>22 . Other authors use the long-time 
self-diffusion constant of a colloid in the dense suspension .D^-iiLIS And in some studies 
the short-time diffusion constant D$ is used, which refers to time-scales, on which a par- 
ticle interacts hydrodynamically but not directly with its effectively static neighbours'^. 
In addition, often D$ and Dl are not measured directly, but inferred from Dq by means 
of two empirical relations: D^/Dq = (1 — <f)/<pg) 2 ' 57 , where <p is the volume fraction of col- 
loidal spheres in the system and 9 is their volume fraction at the glass transition*^ 1 ^, 
and D s /D = (1 — 0/0 rcp ) 117 , where rcp is the volume fraction of spheres at random close 
packing^. We will present our nucleation rate data in terms of all three time-scales for 
comparison, and we will conclude that none of these time-scales suffices to factor out solvent 
effects. 

We have simulated 8240 hard spheres suspended in a solvent at a volume fraction of = 
0.539. The solvent was modelled by means of multi-particle collision-dynamics^— (see 
supplementary material for details). The starting configurations were prepared in the super- 
saturated liquid state and we verified that they did not contain crystalline precursors. Then 
we simulated 20-30 independent trajectories per value of viscosity until crystallization was 
reached in all cases. Despite the relatively high packing fraction, we observe an induction 
period that is long compared to the diffusion times of the system followed by a regime of 
rapid growth (see crystallization curves in the supplementary material). We did not use any 
rare event sampling scheme to speed up sampling, as we did not wish to make any assump- 
tions on the evolution of the density of states or the length of correlation-times involved in 
the process. 

We present all data in units of the sphere diameter a, the mass of a sphere m and the thermal 
energy fcgT. Solvent viscosities ranged between 4 ^fink^T /a 2 and 70 y/mk^T '/a 2 . Translated 
to an experimental system with colloidal particles of radius ~ 420 nm suspended in a solvent 
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FIG. 1. Nucleation rates as a function of solvent viscosity. 

of mass density ~ 1 kg m~ 3 at room temperature, these viscosities correspond to a range of 
8.9 • 10 -6 Pa ■ s to 1.5 • 10~ 4 Pa ■ s. For higher viscosities (up to 2 ■ 10~ 3 Pa ■ s) we simulated the 
diffusion behaviour, but not the crystallization process because the induction times became 
forbiddingling long. 25 

Fig. [1] shows the nucleation rate versus viscosity. The rate drops with increasing viscosity, 
because the colloids' diffusive motion slows down. The question is, whether it is possible to 
rescale the data with respect to the diffusion time such that it collapses onto one value. 

Fig. [2] shows the dependence of the diffusion times on viscosity. We obtained D Q by means 
of a fit to the mean squared displacement of a single particle in the infinitely dilute solvent. 
To determine Dl and D$ we recorded the mean squared displacement of a particle in the 
supersaturated melt. There was a clear diffusive regime at long times to which Dl could 
be fitted. However, even at viscosities that correspond to the range that is commonly 
used in experiments, there was no short-time diffusive regime (for details see supplementary 
material). We therefore resorted to using the maximum value of the time-integral over the 
velocity- velocity autocorrelation function 2 ^ and called it a "short-time diffusion constant", 
although, strictly speaking, it does not correspond to diffusive motion (see supplementary 
material) . 
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FIG. 2. Diffusion constants versus viscosity. Left: Absolute values for Dq (circles), D$ (squares) 

and Dl (stars), Right: Ratios Ds/Dq (triangles), Dl/Dq (squares) and expressions according to 
refi&21. 



The left panel of Fig. [2] shows the absolute values of the diffusion constants, the right panel 
their ratios in comparison to the power laws that are commonly used to estimate them. 
Clearly, the system is not concentrated enough for these assymptotic expressions to hold. 
Nevertheless they have been applied at even lower packing fractions in various studies on 
crystal nucleation in the paslr^ 1 ^. 

Fig. [3] shows the nucleation rates rescaled with respect to Dl and Dq. They clearly depend 
on viscosity. Thus nucleation rates measured in suspensions of hard spheres with different 
solvent viscosities cannot be simply superposed by means of a rescaling with respect to D or 
Dl; crystallization kinetics is clearly influenced by solvent hydrodynamics beyond the level 
of diffusion. 

The lower panel of fig. H] shows the nucleation rates re-scaled with respect to D s . At first 
glance, one might conclude, that this time-scale is the proper choice. However, fig. H] upper 
panel shows a test of this hypothesis for another value of supersaturation (0 = 0.544.). The 
collapse of the data at <fi = 0.539 is a pure coincidence. 

If the nucleation rates depend on viscosity, one could expect to observe differences in the 
sizes, shapes or structures of the crystallites that form, too. We analysed the structures of the 
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FIG. 3. Nucleation rates scaled by Dg (left panel) and Z?l (right panel) as a function of solvent 
viscosity. Both graphs show that the nucleation kinetics change with viscosity beyond the trivial 
dependence of diffusion times. 
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FIG. 4. Nucleation rate normalized with respect to the short-time diffusion constant D s for 
0.539 (lower panel) and (ft = 0.544 (upper panel) 



growing crystallites in terms of their q6q6-bond-order— >2&. Comparing clusters of equal sizes 
from trajectories at different viscosities, neither the radii of gyration nor the distributions of 
q6q6-values of differed within the statistical accuracy. We performed a committor analysis 
for the highest and the lowest viscosity and did not find any difference in the critical cluster 
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size, shape or structure, either. Thus we conclude that the effect we observe is purely kinetic. 

In summary, we have simulated crystallization of a supersaturated suspension of hard spheres 
fully taking into account the solvent hydrodynamics. We find that kinetics need to be taken 
with care when one studies phase transitions in colloids. The crystallization process in hard 
sphere suspensions is richer than has been suggested in the literature so far. There is a clear 
dependence of the nucleation rates on viscosity beyond the trivial slowing down of diffusion. 
However, the crystal structure and shape of the critical nucleus seem to be independent of 
viscosity, hence what we have described here is a purely dynamical effect. In addition we 
showed that there is no clear short-time diffusive regime for the concentrations and viscosities 
commonly used in experiments. And that the power laws that are often used to infer diffusion 
constants do not hold for the volume fraction at which one usually studies nucleation. 

METHODS 

SIMULATION METHOD 

To simulate hard spheres suspeded in a liquid, we used a combination of an event-driven 
molecular dynamics (EDMD) algorithm^ - — for the spheres and multiparticle-collision dy- 
namics (MPCD)22>22 as a mesoscopic solvent model to account for the hydrodynamic interac- 
tions. The basic idea of a MPCD algorithm is to transport momentum through the system 
by means of n point particles of mass m while satisfying the conservation laws of mass, 
energy and momentum locally. The algorithm consists of two steps, namely free streaming 
interrupted by multiparticle collisions. In the streaming step, all fluid particles are propa- 
gated ballistically with their velocities Vjt, i.e. the particle positions are updated according 
to a time increment h as 



After the time step h, the n fluid particles are sorted into a lattice of cubic cells of size 
a x a x a, such that on average n particles are within each collision cell. Then, the particle 



Xfc(t + h) = Xfc(i) + hv k (t). 
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velocities are rotated around the center of mass velocity v in this cell, 

v fc (t + h)= v(t) + 0(a) [v fc (t) - v(*)] , (2) 

where is a rotation matrix corresponding to a fixed angle a, which is generated ran- 

domly for each cell. Before a collision step is carried out, the collision cell grid is shifted by 
a randomly chosen vector with components taken from the interval [—a/2, a/2], to ensure 
Galilean invar iance^. Here, we realized coupling between the colloidal particles and the 
"solvent" particles in the simplest manner, i. e. the colloid takes part in the collision step 
eq. [2] as a point particle with instantaneous velocity V and mass M within its cell with 
v = (MV + m YT v) / (M + nm). 

In order to measure the solvent viscosities we imposed a Poiseuille flow between two planar 
walls. From the resulting parabolic velocity field we extracted 770 from f max - 



SIMULATION DETAILS 

We simulated 8240 hard spheres suspended in a solvent at volume fractions of = 0.539 and 
(ft = 0.544. The starting configurations were prepared in the supersaturated liquid state and 
we verified by means of bond-order analysis that they did not contain cyrstalline precursors. 
We simulated 20-30 independent trajectories per value of viscosity until crystallization was 
reached in all cases. The simulations were run in an unbiased way, i. e. we did not enhance 
statistics by any type of rare event sampling method. 

ANALYSIS 
Nucleation Time 

We identified crystallites by means of the "q6q6-bond order parameter"— 1 ^. Fig |5] shows the 
number of particles in the largest crystalline cluster versus time for two typical crystallization 
runs - one at a low viscosity and one at a high viscosity. The crystallization process clearly 
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FIG. 5. Number of particles in the largest cluster as a function of time for one trajectory at a low 
value of solvent viscosity and one at a high value. 

consists of a very long induction time followed by a fast event. Hence, even though the 
packing fraction is too high to expect classical nucleation theory to hold, we are still dealing 
with an activated, rare event. 

Once a run had produced a cluster of more than 80 particles, it definitely crystallized. Thus 
we used this value to locate the nucleation time. To test the validity of this criterion, 
we performed a commitor analysis for two values of viscosity (77 = 4.17(mfcBT) 1/,2 /a 2 and 
r] = 63.93(m/c B T) 1 / 2 /a 2 ). We found that a cluster size of ca. 30 particles corresponds to a 
50% probability for subsequent full crystallization in both cases. 

As the average induction time we took the arithmetic mean of the distribution of measured 
induction times, and its standard deviation to determine the error bars. The nucleation rate 
is then given by 

1 = V(t ind ) ' 

where V is the volume of the system. 
Diffusion constants 

To determine Do we computed the mean squared displacement of a particle in the infinitely 
dilute solution versus time and fitted a straight line to it. (The result was consistent with 
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FIG. 6. Mean squared displacement of a particle in the suspension versus time for r/ = 
A.n{mkBT) 1 / 2 / a 2 and r/ = 63.93(mA;Br) 1 / 2 /a 2 (main panel) and 915 \fmkoT /a 2 (inset). Even 
at 915 /a 2 which corresponds to 2 • 10~ 3 Pa ■ s there is no short-time diffusive regime. 

the shear viscosity r] that we had measured independently, see above. I.e. the Stokes- 
Einstein relation was fullfilled.) was determined by means of a fit to the mean squared 
displacement of a particle in the dense solution (see fig. E]). Even for the highest viscosity 
we simulated, which corresponds to the viscosity of a typical experiment on hard sphere 
crystallization, we did not see a short time diffusive regime. 

In order to obtain a value for D$, nevertheless, we tried two other approaches: we fitted the 
dynamic structure factor, but the results dependend strongly on the choice of time-window 
for the "short-time" regime. And we computed the time-integral over the velocity-velocity 
autocorrelation function, as suggested in ref.— 



and identified its maximum value with D$ (see Fig. [7j). This approach turned out to be 
relativley robust. 

Structure of Crystallites 
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FIG. 7. Measured time integrals over velocity velocity autocorrelation function. The maxima were 
taken for Dg 
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FIG. 8. Distributions of q^qe for different sizes of the largest cluster. Here shown for small and 
high viscosity. 



To check if the structure of the crystallites depends on viscosity, we computed the qeq^ dis- 
tributions for clusters of equal sizes at different viscosities. Fig. |H] shows that the crystallites 
are very similar in structure. The snapshots figs. [91 [10]do not allow for a distinction between 
low and high viscosity, either. 
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FIG. 9. Timeseries of snapshots of the largest cluster for small viscosity. 




FIG. 10. Timeseries of snapshots of the largest cluster for high viscosity. 
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